home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d18 / nrpas13.arc / ANNEAL.PAS < prev    next >
Pascal/Delphi Source File  |  1991-05-01  |  5KB  |  172 lines
  1. PROCEDURE anneal(x,y : cityarray; VAR iorder: iarray; ncity: integer);
  2. (* Programs using routine ANNEAL must define types
  3.    cityarray : ARRAY [1..ncity] OF real;
  4.    iarray : ARRAY [1..ncity] OF integer;
  5. in the main routine. *)
  6. LABEL 10,20,99;
  7. CONST
  8.    tfactr = 0.9;
  9. TYPE
  10.    nsix = ARRAY [1..6] OF integer;
  11. VAR
  12.    ans : boolean;
  13.    path,de,t : real;
  14.    nover,nlimit,i1,i2,idum,iseed: integer;
  15.    i,j,k,nsucc,nn,idec : integer;
  16.    n : nsix;
  17.  
  18. FUNCTION alen(x1,x2,y1,y2: real): real;
  19. BEGIN
  20.    alen := sqrt(sqr(x2-x1)+sqr(y2-y1))
  21. END;
  22.  
  23. PROCEDURE revcst(x,y: cityarray; iorder: iarray; ncity: integer;
  24.       VAR n: nsix; VAR de: real);
  25. VAR
  26.    xx,yy : ARRAY [1..6] OF real;
  27.    j,ii : integer;
  28. BEGIN
  29.    n[3] := 1 + ((n[1]+ncity-2) MOD ncity);
  30.    n[4] := 1 + (n[2] MOD ncity);
  31.    FOR j := 1 TO 4 DO BEGIN
  32.       ii := iorder[n[j]];
  33.       xx[j] := x[ii];
  34.       yy[j] := y[ii]
  35.    END;
  36.    de := -alen(xx[1],xx[3],yy[1],yy[3])-alen(xx[2],xx[4],yy[2],yy[4])
  37.       +alen(xx[1],xx[4],yy[1],yy[4])+alen(xx[2],xx[3],yy[2],yy[3])
  38. END;
  39.  
  40. PROCEDURE reverse(VAR iorder: iarray; ncity: integer; n: nsix);
  41. VAR
  42.    nn,j,k,l,itmp : integer;
  43. BEGIN
  44.    nn := (1+((n[2]-n[1]+ncity) MOD ncity)) DIV 2;
  45.    FOR j := 1 TO nn DO BEGIN
  46.       k := 1 + ((n[1]+j-2) MOD ncity);
  47.       l := 1 + ((n[2]-j+ncity) MOD ncity);
  48.       itmp := iorder[k];
  49.       iorder[k] := iorder[l];
  50.       iorder[l] := itmp
  51.    END
  52. END;
  53.  
  54. PROCEDURE trncst(x,y: cityarray; iorder: iarray; ncity: integer;
  55.       VAR n: nsix; VAR de: real);
  56. VAR
  57.    xx,yy : ARRAY [1..6] OF real;
  58.    j,ii : integer;
  59. BEGIN
  60.    n[4] := 1 + (n[3] MOD ncity);
  61.    n[5] := 1 + ((n[1]+ncity-2) MOD ncity);
  62.    n[6] := 1 + (n[2] MOD ncity);
  63.    FOR j := 1 TO 6 DO BEGIN
  64.       ii := iorder[n[j]];
  65.       xx[j] := x[ii];
  66.       yy[j] := y[ii]
  67.    END;
  68.    de := -alen(xx[2],xx[6],yy[2],yy[6])-alen(xx[1],xx[5],yy[1],yy[5])
  69.       -alen(xx[3],xx[4],yy[3],yy[4])+alen(xx[1],xx[3],yy[1],yy[3])
  70.       +alen(xx[2],xx[4],yy[2],yy[4])+alen(xx[5],xx[6],yy[5],yy[6])
  71. END;
  72.  
  73. PROCEDURE trnspt(VAR iorder: iarray; ncity: integer; n: nsix);
  74. CONST
  75.    maxcity=1000;
  76. VAR
  77.    jorder : ARRAY [1..maxcity] OF integer;
  78.    m1,m2,m3,nn,j,jj : integer;
  79. BEGIN
  80.    m1 := 1 + ((n[2]-n[1]+ncity) MOD ncity);
  81.    m2 := 1 + ((n[5]-n[4]+ncity) MOD ncity);
  82.    m3 := 1 + ((n[3]-n[6]+ncity) MOD ncity);
  83.    nn := 1;
  84.    FOR j := 1 TO m1 DO BEGIN
  85.       jj := 1 + ((j+n[1]-2) MOD ncity);
  86.       jorder[nn] := iorder[jj];
  87.       nn := nn+1
  88.    END;
  89.    IF (m2>0) THEN BEGIN
  90.       FOR j := 1 TO m2 DO BEGIN
  91.          jj := 1+((j+n[4]-2) MOD ncity);
  92.          jorder[nn] := iorder[jj];
  93.          nn := nn+1
  94.       END
  95.    END;
  96.    IF (m3>0) THEN BEGIN
  97.       FOR j := 1 TO m3 DO BEGIN
  98.          jj := 1 + ((j+n[6]-2) MOD ncity);
  99.          jorder[nn] := iorder[jj];
  100.          nn := nn+1
  101.       END
  102.    END;
  103.    FOR j := 1 TO ncity DO BEGIN
  104.       iorder[j] := jorder[j] 
  105.    END
  106. END;
  107.  
  108. PROCEDURE metrop(de,t: real; VAR ans: boolean);
  109. (* Programs using routine METROP must declare the variable 
  110. VAR
  111.    gljdum : integer;
  112. and initialize its value to
  113.    gljdum := 1;
  114. in the main routine. *)
  115. BEGIN
  116.    ans := (de<0.0) OR (ran3(gljdum)<exp(-de/t))
  117. END;
  118.  
  119. BEGIN
  120.    nover := 100*ncity;
  121.    nlimit := 10*ncity;
  122.    path := 0.0;
  123.    t := 0.5;
  124.    FOR i := 1 TO (ncity-1) DO BEGIN
  125.       i1 := iorder[i];
  126.       i2 := iorder[i+1];
  127.       path := path+alen(x[i1],x[i2],y[i1],y[i2])
  128.    END;
  129.    i1 := iorder[ncity];
  130.    i2 := iorder[1];
  131.    path := path+alen(x[i1],x[i2],y[i1],y[i2]);
  132.    idum := -1;
  133.    iseed := 111;
  134.    FOR j := 1 TO 100 DO BEGIN
  135.       nsucc := 0;
  136.       FOR k := 1 TO nover DO BEGIN
  137. 10:      n[1] := 1+trunc(ncity*ran3(idum));
  138.          n[2] := 1+trunc((ncity-1)*ran3(idum));
  139.          IF (n[2]>=n[1]) THEN n[2] := n[2]+1;
  140.          nn := 1+((n[1]-n[2]+ncity-1) MOD ncity);
  141.          IF (nn<3) THEN goto 10;
  142.          idec := irbit1(iseed);
  143.          IF (idec=0) THEN BEGIN
  144.             n[3] := n[2]+trunc(abs(nn-2)*ran3(idum))+1;
  145.             n[3] := 1+((n[3]-1) MOD ncity);
  146.             trncst(x,y,iorder,ncity,n,de);
  147.             metrop(de,t,ans);
  148.             IF ans THEN BEGIN
  149.                nsucc := nsucc+1;
  150.                path := path+de;
  151.                trnspt(iorder,ncity,n)
  152.             END
  153.          END ELSE BEGIN
  154.             revcst(x,y,iorder,ncity,n,de);
  155.             metrop(de,t,ans);
  156.             IF ans THEN BEGIN
  157.                nsucc := nsucc+1;
  158.                path := path+de;
  159.                reverse(iorder,ncity,n)
  160.             END
  161.          END;
  162.          IF (nsucc>=nlimit) THEN goto 20
  163.       END;
  164. 20:   writeln;
  165.       writeln('T =',t:10:6,'    Path Length =',path:12:6);
  166.       writeln('Successful Moves: ',nsucc:6);
  167.       t := t*tfactr;
  168.       IF (nsucc=0) THEN goto 99
  169.    END;
  170. 99:
  171. END;
  172.